
############################################
############################################
# previous: 8-VAP_CDCut ####################
############################################

###### Summary ############################################################################################
#This script merges all species for each 'merging unit' together: A merging unit represents closely related species 
#with similar venom composition and similar maximum abundances, which need to be treated as one continuum of habitat 
#suitability rather than separate, overlapping taxonomic units, ie. they are not additive to each other but rather 
#'flow' from one species to the next and represent a continuous distribution that affects victims in a similar fashion.
#merging maps are created by calculating the maximum habitat suitability across all models in the merging groups and
#show groups of species that...
#are closely related AND have similar venom properties AND ...
#are likely to have one dominant species per pixel,i.e. range overlaps don't mean that they actually co-occur but rather 
#represent uncertainty of which species is present in transition zones OR for which range overlap might be true but one 
#species is likely to greatly outcompete the others in each pixel, making the main species' possible maximum abundance 
#representative of total abundance for thatgroup in that pixel.
#Once merging units are created, the second part of the script then adds their binary or continuous distributions together 
#to create hotspot maps: maps that show diversity and cumulative abundance for different 'hotspot gorups' such as all snakes, all 
#vipers,all merging groups that contain category one species, all snakes with a certain venom activity, all snakes covered by as 
#certain polyvalent antivenom, snakes that require treatment within x hours, etc. The results are plotted as .png images.
###########################################################################################################

############################################
# next: 10_VAP_ClimatChange_Summaries ######
############################################
#other potential next scripts: 
#create per administrative unit summaries of outputs to use in snakebite incidence models
#weight hotspots by actual abundance (suitability only shows abundance relative to maximum abundance 0-1), aggressiveness, etc.
#calculate likely snakebite incidence from first principals based on abundance, human exposure, etc.
############################################

#################################################################################################################
#################################################################################################################
#################################################################################################################

#Packages:
library(R.utils)
library(raster)
library(maptools)
library(rgdal)
library(stringr)
library(sf)

#Project Info:
#Project name:
projectName<-"WHOSnakes"
#Resolution in decimal degrees:
res<- 0.01
#define selection criteria of interest for further analyses:
crit<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")
critThresh<-c(1)                     #c(1,1,1,10,1,75)

#Directories:
MyDir<-paste("/home/jc217070/Maxent",sep="")
OccDir<-file.path(MyDir,"Maxent_SWDs/spp_occs")
FinalOutDir<-paste(MyDir,"/Outputs3_Final",sep="")
HotSpotDir<-paste(FinalOutDir, "/HotSpots",sep="")
IndSppOutDir<-paste(FinalOutDir, "/Final_Individual_Spp",sep="")
MergeSppOutDir<-paste(FinalOutDir, "/Final_Merge_Groups",sep="")

PlotDir<-paste(MyDir,"/Maxent_Data/Plots/Plots_WHOSnakes_2020_09_01",sep="")

#Species Info:
#DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,"_Fin3.csv",sep="")), header=TRUE, sep=",") #load data summary for all species
spp<-DatSum$gen_sp_subtax

#Spatial Data:
worldmap <- st_read(paste(PlotDir,"/worldmap/WHO_countries.shp",sep=""))
WHOpolies <- st_read(paste(MyDir,"/WHOSnakesShapes/SSN_CTY_CAT_BRG_FEA_IN_Mar2023.shp",sep=""))

bg<-raster(paste(PlotDir,"/Cop_FAPAR_mean_2010-2019_WW1km.asc",sep=""), crs="+init=epsg:4326")
GLOBraster<-bg
values(GLOBraster)<-ifelse(values(GLOBraster)!=-Inf,0,-Inf) #create 'empty' raster of same extent/res as FAPAR
GLOBrasterB<-GLOBraster
values(GLOBrasterB)<-ifelse(values(GLOBrasterB)==0,1,values(GLOBrasterB)) #create 'empty' raster of same extent/res as hillshade

my_col1 = rev(heat.colors(n = 100, alpha=1))
my_col2<-rgb(74,122,218, max = 255, alpha = 255)
my_col3<-rgb(74,122,218, max = 255, alpha = 100)
mycols1<-c("white","grey")
mycols2<-c("royalblue","cyan","navy","skyblue","black","gold",
           "red","darkred","darkorange","orangered","deeppink","pink",
           "purple","slateblue1","lightsalmon","burlywood4","coral3",
           "deeppink4","darkorange4","cyan4","purple4")
WHOcol <- rgb(150, 180, 250, max = 255, alpha = 0) #alpha=0 means completely transparent

#####################################################################################
#####################################################################################
#####################################################################################
#####################################################################################
#####################################################################################


# B.
#####################################################################################
# make hotspot maps of species diveristy and cumulative suitability for this group ##
# of species (e.g. all SEA spitting Cobras) for each current & future category ######
#####################################################################################
########################
#################################################################
#define names of groups to make hotspots for using merged rasters:
#e.g. groupnamesMerge<-c("all", "includescat1","vipers","elapids","colubrids","includescat1Nigeria")
#e.g. groupnamesMergeFull<-c("All medically relevant snakes (grouped)","all cat 1 snakes (grouped)",
#                       "all medicaly relevant vipers (grouped)","all medically relevant elapids (grouped)",
#                       "all medically relevant colubrids (grouped)","groups present in Nigeria including category 1 snakes (grouped)")
groupnamesMerge<-c("all")
groupnamesMergeFull<-c("All medically relevant snakes (grouped)")

#define species groups included in each of these groups:
# includescat1Nigeria<-unique(DatSum$MergeShort[DatSum$RiskCat==1 & DatSum$gen_sp_subtax %in% (DatSum$gen_sp_subtax[grep("NGA", DatSum$country, fixed = TRUE)])])
# selectNigeria<-c("Egyptian Cobra Relatives (Uraeus complex)","African True Spitting Cobras (Afronaja complex)",
#                  "Forest cobras (Boulengerina melanoleuca group)","Green Mambas (Dendroaspis viridis group)",
#                  "Black Mamba (Dendroaspis polylepis)","Boomslang (Dispholidus typus)",
#                  "Western Carpet Vipers (Echis coloratus & ocellatus complex)","Horned Vipers (Cerastes spp.)",
#                  "Puff Adder (Bitis arietans)")
# allNigeria<-unique(DatSum$MergeShort[DatSum$RiskCat==1 & DatSum$gen_sp_subtax %in% (DatSum$gen_sp_subtax[grep("NGA", DatSum$country, fixed = TRUE)])])
all<-unique(DatSum$MergeShort[DatSum$Listing_comment!="not listed"])
  
#list species lists for all these groups:
#groupnamesMerge contains the group names and 
#groupsppMerge contains the associated species lists in the same order
groupsppMergeRaw<-list(all)
groupsppMerge<-list(gsub(" ","_",gsub("\\.","",all)))
###################################################################

#################################################################
#define names of groups to make hotspots for using individual species rasters:
# groupnamesIndiv<-c("Big4India","onlycat1Nigeria","onlycat2Nigeria")
# groupnamesIndivFull<-c("Big 4 of India (individual species)","category 1 snakes present in Nigeria (individual species)",
#                        "category 2 snakes present in Nigeria (individual species")
groupnamesIndiv<-c("Big4India")
groupnamesIndivFull<-c("Big 4 of India (individual species)")

#define species included in each of these groups:
#groupnamesIndiv contains the group names and 
#groupsppIndiv contains the associated species lists in the same order
Big4India<-c("Bungarus_caeruleus","Echis_carinatus","Naja_naja","Daboia_russelii")
#onlycat1Nigeria<-DatSum$gen_sp_subtax[DatSum$RiskCat==1 & DatSum$gen_sp_subtax %in% (DatSum$gen_sp_subtax[grep("NGA", DatSum$country, fixed = TRUE)])]
#onlycat2Nigeria<-DatSum$gen_sp_subtax[DatSum$RiskCat==2 & DatSum$gen_sp_subtax %in% (DatSum$gen_sp_subtax[grep("NGA", DatSum$country, fixed = TRUE)])]
#make list object of species groups for all these hotspot groups:
#groupsppIndiv<-list(Big4India,onlycat1Nigeria,onlycat2Nigeria)
groupsppIndiv<-list(Big4India)

###################################################################

###################################################################
#make final list of all groups, those from merged rasters 
#and those from individual species rasters:
groupnames<-c(groupnamesMerge,groupnamesIndiv)
groupnamesFull<-c(groupnamesMergeFull,groupnamesIndivFull)
###################################################################

HotSpotCats<-c("Suit","Bin")
#HotSpotCats<-c("Bin")

OutCats<-c("Current","2050_Median","2090_Median","2090_Q10","2090_Q90","2050_Q10","2050_Q90")
#OutCats<-c("2050_Median")
########################

Sys.sleep(0.1)
print(paste(Sys.time(), "making hotspot maps of species and/or species groups", sep=" "))

#####################################################################
for(groupname in groupnames[1]){
  
  # #################################################################################
  # #make temporary temp dir for rasters in this job and set from default to this dir
  # dir.create(paste(FinalOutDir,"/",groupname,"_tmpHotSpot",sep=""),showWarnings=FALSE)
  # rasterOptions(tmpdir=paste(FinalOutDir,"/",groupname,"_tmpHotSpot",sep=""))
  # #################################################################################
  # 
  k<-which(groupnames==groupname)
  Sys.sleep(0.1)
  print(paste(Sys.time(), "making merged maps for",groupname, sep=" "))
  
  if(groupname %in% groupnamesMerge){
    thisgrouplist<-"groupnamesMerge"
    i<-which(groupnamesMerge==groupname)
    myspp<-groupsppMerge[[i]]
    RasterDir<-"/home/jc217070/Maxent/Outputs3_Final/Final_Merge_Groups"
  }
  
  if(groupname %in% groupnamesIndiv){
    thisgrouplist<-"groupnamesIndiv"
    i<-which(groupnamesIndiv==groupname)
    print(paste(groupnamesIndiv[[i]]))
    myspp<-groupsppIndiv[[i]]
    RasterDir<-"/home/jc217070/Maxent/Outputs3_Final/Final_Individual_Spp"
  }
  
  for(cat in OutCats){ #for each current & future category
    for(BinOrSuit in HotSpotCats){ #for binary and continuous outputs
      
      removeTmpFiles(h=0.5) #removes all temporary files older than 30min

      #change the beginning species for the last failed 
      #batch to the point where it needs to restart from
      ####################################      
      if(BinOrSuit=="Suit"){
        j<-1
      } else {
        j<-1
      }
      ####################################      

      HotSpot<-GLOBraster
      
      for (sp in myspp [j:length(myspp)]){ #for each species in this subset

        n<-which(myspp==sp)[1]
        removeTmpFiles(h=0.5) #removes all temporary files older than 30min, then read hotspot map back in to continue process
        
        if(thisgrouplist=="groupnamesIndiv"){
          if (BinOrSuit=="Bin"){
            outfile<-paste(RasterDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc",sep="")
          }else{
            outfile<-paste(RasterDir,"/",sp,"_",cat,"_CropThreshCDCut.asc",sep="")
          }
        }
        
        if(thisgrouplist=="groupnamesMerge"){
          outfile<-paste(RasterDir,"/",sp,"_",cat,BinOrSuit,"_Merged.asc",sep="")
        }
        
        if (file.exists(paste(outfile,".gz", sep="")) & !(file.exists(paste(outfile)))){ #safety, make sure that output .asc exists
          
          gunzip(filename=paste(outfile,".gz", sep=""), #unzip model
                 destname=paste(outfile),
                 remove=FALSE, overwrite=TRUE)
        }
        
        if (file.exists(paste(outfile, sep=""))){ #safety, make sure that output .asc exists
          
          #print out summary of whats going on (keep track of progress: start)
          Sys.sleep(0.1)
          print(paste(Sys.time(), "adding", sp, "to merge map, j=",j, sep=" "))
          
          #unzip & read thresholded sdm
          SDM<-raster(paste(outfile), crs="+init=epsg:4326") #read model output
          SDM<-merge(SDM,GLOBraster) #merge to GLOBraster to make all rasters sae extent for adding together
          
          
          ############################
          # add to previous SDMs to create hotspot map
          if(file.exists(paste(HotSpotDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_HotSpot",j-1,".asc",sep=""))){
            HotSpot<-raster(paste(HotSpotDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_HotSpot",j-1,".asc",sep=""), crs="+init=epsg:4326") #read model output for species
          }
          
          HotSpot<-HotSpot+SDM
          writeRaster(HotSpot, filename=paste(HotSpotDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_HotSpot",j,".asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
          
          #remove previous backup file
          if(file.exists(paste(HotSpotDir,"/",groupname,"_",cat, BinOrSuit,"_prelim_HotSpot",j-2,".asc",sep=""))){
            file.remove(paste(HotSpotDir,"/",groupname,"_",cat, BinOrSuit,"_prelim_HotSpot",j-2,".asc",sep=""))
          }

          if (!file.exists(paste(outfile,".gz", sep="")) & (file.exists(paste(outfile)))){ #safety, make sure that output .asc exists
            gzip(filename=paste(outfile),overwrite=TRUE, remove=TRUE)
          }
          
          if (file.exists(paste(outfile,".gz", sep="")) & (file.exists(paste(outfile)))){ #safety, make sure that output .asc exists
            file.remove(filename=paste(outfile))
          }
          ############################
          
          
          
        }  else { #end safety loop
          
          #print out summary of whats going on (keep track of progress: start)
          Sys.sleep(0.1)
          print(paste(Sys.time(), " can't find final output for ", sp, ", not added to hotspot map", sep=""))
        }
        
        j<-j+1
        
      } #end species loop
      
      Sys.sleep(0.1)
      print(paste(Sys.time(), " done with species loop for ", groupname, " ", BinOrSuit, sep=""))
      
      
      if(file.exists(paste(HotSpotDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_HotSpot",j-1,".asc",sep=""))){
        HotSpot<-raster(paste(HotSpotDir,"/",groupname,"_",cat,BinOrSuit,"_prelim_HotSpot",j-1,".asc",sep="")) #read model output for species
      }
      
      
      #remove previous backup file
      if(file.exists(paste(HotSpotDir,"/",groupname,"_",cat, BinOrSuit,"_prelim_HotSpot",j-2,".asc",sep=""))){
        file.remove(paste(HotSpotDir,"/",groupname,"_",cat, BinOrSuit,"_prelim_HotSpot",j-2,".asc",sep=""))
      }
      
      
      maprecs<-data.frame()
      
      if(groupname %in% groupnamesMerge){
        i<-which(groupnamesMerge==groupname)
        for (sp in myspp){ #for each species in this subset
          for (taxon in DatSum$gen_sp_subtax[DatSum$MergeShort %in% groupsppMergeRaw[[i]]]){
            if(file.exists(paste(OccDir,"/",taxon,"_allRecs.csv",sep=""))){
              #get species records used for model:
            addrecs<-read.table(paste(OccDir,"/",taxon,"_allRecs.csv",sep=""), sep=",", header=TRUE)   #read in species recs
            maprecs<-rbind(maprecs,addrecs)
          }
        }
      }}
      
      if(groupname %in% groupnamesIndiv){
        for (sp in myspp){ #for each species in this subset
          if(file.exists(paste(OccDir,"/",sp,"_allRecs.csv",sep=""))){
            #get species records used for model:
          addrecs<-read.table(paste(OccDir,"/",sp,"_allRecs.csv",sep=""), sep=",", header=TRUE)   #read in species recs
          maprecs<-rbind(maprecs,addrecs)
        }
      }}
      
      
      if(length(maprecs[,1])>0){
        
        Sys.sleep(0.1)
        print(paste(Sys.time(), "making extent for modelling unit ",groupname, sep=" "))
        
        rangewidthDG<-max((max(na.omit(maprecs$long))-min(na.omit(maprecs$long))),
                          (max(na.omit(maprecs$lat))-min(na.omit(maprecs$lat))))
        rangewidthM<-rangewidthDG*100000 #turn dec deg into meters
        rangewidthM<-ifelse(rangewidthM>3000000,3000000,rangewidthM) #make buffer no bigger than this
        rangewidthM<-ifelse(rangewidthM<1000000,1000000,rangewidthM) #make buffer no smaller than this
        y<-rangewidthM

        #make extent object describing a box around occurrence records (margin= 4 times x and y extent of occurrence records):
        xrange<- (max(na.omit(maprecs$long)))-(min(na.omit(maprecs$long)))
        yrange<- (max(na.omit(maprecs$lat)))-(min(na.omit(maprecs$lat)))
        xcenter<- (min(na.omit(maprecs$long)))+(xrange/2)
        ycenter<- (min(na.omit(maprecs$lat)))+(yrange/2)
        rangemax<-max(xrange,yrange)
        xmin<-ifelse((xcenter-(rangemax/2)- y/100000) >= -180, (xcenter-(rangemax/2)- y/100000), -180)
        xmax<-ifelse((xcenter+(rangemax/2)+ y/100000) <=  180, (xcenter+(rangemax/2)+ y/100000),  180)
        ymin<-ifelse((ycenter-(rangemax/2)- y/100000) >=  -90, (ycenter-(rangemax/2)- y/100000),  -90)
        ymax<-ifelse((ycenter+(rangemax/2)+ y/100000) <=   90, (ycenter+(rangemax/2)+ y/100000),   90)
        myext<-extent(xmin,xmax,ymin, ymax)
      } else{
        myext<-extent(-180,180,-90,90)
      }
      
      HotSpot<-crop(HotSpot,myext)
      
      #write final version:
      writeRaster(HotSpot, filename=paste(HotSpotDir,"/",groupname,"_",cat,BinOrSuit,"_HotSpot.asc",sep=""), format="ascii", overwrite=TRUE) #write thresholded model raster
      removeTmpFiles(h=0.5)
      
      
      ########################################################################
      #remove any temporary merge-files after making all hotspot maps:
      temp.file.list<-list.files(path=paste(HotSpotDir),pattern=paste(groupname,"_",cat,BinOrSuit,"_prelim_HotSpot",sep=""),full.names=TRUE)
      for (file in temp.file.list){
        file.remove(paste(file))
      }
    }
    #######################################################################
    
    ################################################################################
    ### plot maps for this subset ##################################################
    
    Sys.sleep(0.1)
    print(paste(Sys.time(), ", plotting hotspot map", sep=" "))
    
    HotSpotBin<-raster(x=paste(HotSpotDir,"/",groupname, "_",cat,"Bin_HotSpot.asc",sep=""), crs="+init=epsg:4326") #read subset hotspot back in
    HotSpotSuit<-raster(x=paste(HotSpotDir,"/",groupname, "_",cat,"Suit_HotSpot.asc",sep=""), crs="+init=epsg:4326") #read subset hotspot back in
    
    png(filename=paste(PlotDir,"/",groupname,"_",cat,"_HotSpot.png",sep=""), units="in", width=10, height=14, res=500) #png smaller file size, A3 for better resolution
    par(mfrow = c(2,1), mar = c(1,1,0,0), oma=c(1,1,4,1), mgp=c(1.2,0.5,0))
    
    plot(HotSpotSuit, col=my_col1, box=FALSE, bty="n", axes=FALSE, legend=TRUE)
    plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
    
    plot(HotSpotBin, col=my_col1, box=FALSE, bty="n", axes=FALSE, legend=TRUE)
    plot(worldmap[5], border="grey20", lwd= 0.5, add=TRUE,col=NA,reset=FALSE)
    
    mtext(paste(cat," cumulative relative abundance (using suitability as a proxy; top)",sep=""),
          side = 3, line = 2, outer = TRUE, cex= 1.3)
    mtext(paste("& number of species present (bottom) for",sep=""), 
          side = 3, line = 1, outer = TRUE, cex= 1.3)
    mtext(paste("'",groupnamesFull[k],"'",sep=""), 
          side = 3, line = 0, outer = TRUE, cex= 1.3)
    
    
    dev.off()
    
    Sys.sleep(0.1)
    print(paste(Sys.time(), ", done plotting hotspot map", sep=" "))
    ########################################################################
   
    #zip final files
    if(file.exists (paste(HotSpotDir,"/",groupname,"_",cat,"Bin_HotSpot.asc",sep=""))){
      gzip(filename=paste(HotSpotDir,"/",groupname,"_",cat,"Bin_HotSpot.asc",sep=""),overwrite=TRUE, remove=TRUE)
    }
    
    if(file.exists (paste(HotSpotDir,"/",groupname,"_",cat,"Suit_HotSpot.asc",sep=""))){
      gzip(filename=paste(HotSpotDir,"/",groupname,"_",cat,"Suit_HotSpot.asc",sep=""),overwrite=TRUE, remove=TRUE)
    }
    
     
  }
  
  # ######################################################
  # #delete the temp dir for this script
  # unlink(paste(FinalOutDir,"/",groupname,"_tmpHotSpot",sep=""), recursive = TRUE) 
  # ######################################################
}

################################################################################
################################################################################
################################################################################
################################################################################
################################################################################
################################################################################